      program tilehamfit
      implicit double precision(a-h,o-z)
      include 'common-thfit.f'
      parameter(lw=2*(7*nn+nn*nn+2*mm*nn+3*mm+nn*(nn-1)),liw=lw)
      parameter(ljc=mm)
c      dimension xc(nn),fvecc(mm)  ! for user-def LSFUN
      dimension a(nn),iw(liw),w(lw) ! for e04fdf
      integer po(2,mm)
      integer nel(nn),dn2(nn),cpo(nn)
      character*50 oname(nn),tname(mm)
      integer iusr(lw),nat(mm)
      double precision usr(lw)
      external LSFUN2
c      call x04aaf(1,7)
      open(2,file='thfit.inp')
      read(2,*)nncols,neli,(nel(io),io=1,neli) !columns, which eliminate
      read(2,*) !(oname(io),io=1,nncols)
      do io=1,nncols
         oname(io)='x'
      enddo
      ndata=0
      do io=1,nn
         dn2(io)=0
      enddo
      do io=1,mm
         read(2,*,end=888)
     +        (xob_(jo,io),jo=1,nncols),ee(io),nat(io),tname(io)
         ndata=ndata+1
         do jo=1,nncols
            dn2(jo)=dn2(jo)+int(xob_(jo,io))**2
         enddo
c         po(1,ndata)=i1
c         po(2,ndata)=i2
c         write(*,*)x,x,(xob(jo,io),jo=1,nncols),ee(io)
      enddo
 888  close(2)
      do io=1,neli
         dn2(nel(io))=0
      enddo
      ncols=0
      do jo=1,nncols ! eliminate unused parameters
         if (dn2(jo).gt.0) then
            ncols=ncols+1
            cpo(ncols)=jo
            do io=1,ndata
               xob(ncols,io)=xob_(jo,io)
            enddo
         endif
      enddo
      do io=1,ncols
         a(io)=0d0
      enddo
c      a(ncols)=-100000d0
      fsumsq=0d0
ccc   NOTE: e04gef was in NAG 87 release. e04gzf is used in NAG 97 release.
c      call e04gef(ndata,ncols,a,fsumsq,iw,liw,w,lw,ifail) !calls LSFUN2
c      a(1)=1d0
c      a(2)=-1d0
c      a(3)=-1d0
c      a(4)=1d0
      call e04gzf(ndata,ncols,LSFUN2,a,fsumsq,w,lw,iusr,usr,ifail) !calls LSFUN2
c      write(*,*)'IFAIL ',ifail
      write(*,223)ndata
      write(*,224)ncols
      write(*,225)(df2/ndata)**0.5
c      open(2,file='objects.out')
c      do io=1,nncols
c         read(2,*)oname(io)
c      enddo
c      close(2)
c      write(*,221)
      do io=1,ncols
         write(*,222)cpo(io),a(io),oname(cpo(io))
      enddo
 221  format(" #  COEF[eV]  OBJECT")
 222  format('coef ',i2,f10.3,'  ',a30)
      open(2,file='thfit.out')
      x2=0
      do io=1,ndata
         cc=0d0
         do jo=1,ncols
            cc=cc+xob(jo,io)*a(jo)
         enddo
         x=ee(io)-cc
         xn=x/nat(io)
         x2=x2+xn*xn
         write(2,227)ee(io),cc,x,xn,nat(io),tname(io) !,po(1,io),po(2,io)
      enddo
      write(*,226)(x2/ndata)**0.5*1000
      close(2)
 224  format('parameters ',i7)
 223  format('datapoints ',i7)
 225  format('rms_scell ',f7.4,' meV')
 226  format('rms_atom ',f10.2,' meV')
 227  format(4f13.5,i6,' ',a20)
      end
c      subroutine LSFUN2(M,N,XC,FVECC,FJACC,LJC)
      subroutine LSFUN2(M,N,XC,FVECC,FJACC,LJC,iusr,usr)
      implicit double precision(a-h,o-z)
      include 'common-thfit.f'
      double precision XC(N) !variables: hamiltonian
      double precision FVECC(M) !function values
      double precision FJACC(LJC,N) !partial derivatives
      integer iusr(*)
      double precision usr(*)
      df2=0d0
      do io=1,m
         calc=0d0
         do jo=1,n
            calc=calc+xob(jo,io)*xc(jo)
            fjacc(io,jo)=xob(jo,io)
         enddo
         df=calc-ee(io)
         fvecc(io)=df
         df2=df2+df*df
      enddo
c      write(*,333)(df2/m)**0.5,(xc(io),io=1,n)
 333  format('End LSF : ',f8.4,' pars ',15f9.5)
      return
      end

